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ABSTRACT 

A fundamental assumption in our understanding of disks is that when the Toomre parameter Q ^> 1, 
£f~) the disk is stable against fragmentation into smaller, self-gravitating objects (and so cannot, for example, 

t— I form planets via direct collapse). However, if disks are turbulent, this criterion neglects a broad spectrum 

of stochastic density fluctuations that can produce rare but extremely high-density local mass concentra- 
tions that will easily collapse. We have recently developed an analytic framework to predict the statistics 
of these fluctuations. Here, we use these models to consider the rate of fragmentation and mass spectrum 
of fragments formed in a turbulent, Keplerian disk (e.g. a proto-planetary or proto-stellar disk). Turbulent 
disks are never completely stable: we calculate the (always finite) probability of the formation of self- 
i— ( gravitating structures via stochastic turbulent density fluctuations (compressions, shocks, etc.) in such a 

t— I disk. Modest sub-sonic turbulence above a minimum Mach number > 0. 1 is sufficient to produce a few 

stochastic fragmentation or "direct collapse" events over ^Myr timescales, even if Q 1 and cooling is 
relatively "slow" (f coo i 3> f rbit)- In trans-sonic turbulence (Ai ~ 1) this extends to Q as large as ~ 100. 
We derive the "true" Q criterion needed to suppress such events (over some timescale of interest), which 
scales exponentially with Mach number. We specify this to cases where the turbulence is powered by MRI, 
f~| convection, or density/spiral waves, and derive equivalent criteria in terms of Q and/or the disk cooling 

Qh time. In the latter case, cooling times as long as > 50 SI 1 may be required to completely suppress this 

q channel for collapse. These gravo-turbulent events produce a mass spectrum concentrated near the Toomre 

mass ~ (Q M / M *) 2 M disk (spanning rocky-to-giant planet masses, and increasing with distance from 
the star), with a wider mass spectrum as Ai increases. We apply this to plausible models of proto-planetary 
disks with self-consistent cooling rates and temperatures, and show that, beyond radii ~ 1 — 10 au, no disk 
1—1 temperature can fully suppress stochastic events. For theoretically expected temperature profiles, even a 

minimum mass solar nebula could experience stochastic collapse events, provided a source of turbulence 
with modest sub-sonic Mach numbers. 

Key words: protoplanetary discs — planets and satellites: formation — accretion, accretion disks — 
hydrodynamics — instabilities — turbulence — star formation: general — galaxies: formation 

1 INTRODUCTION proto-planetary disks survive, to account for a significant fraction 

t—( of planets. 

1.1 The Problem: When Do Disks Fragment? Yet there h &m nQ consensus in the literature regarding the 

Fragmentation and collapse of self-gravitating gas in turbulent criteria for "stability" versus fragmentation in any disk, especially 

^ disks is a process central to a wide range of astrophysics, including in nearly- Keplerian proto-planetary disks. The classic Toomre Q 

^> planet, star, supermassive black hole, and galaxy formation. The criterion: 

specific case of a "marginally stable," modestly turbulent Keple- q_ a z K ^ ^ ^ 



PL, 



x 

hi planetary and proto-stellar disks as well as AGN accretion disks. , . , .... „ , .. . 

Cu L,, ui.ii, -ii u where a g is the gas velocity dispersion, k ~ \l the epicychc ire- 



rian disk is particularly important, since this is expected in proto- 7rG£ g . 
planetary and proto-stellar disks as well as AGN accretion disks. 

There has been considerable debate, especially, regarding whether , , ,. , , „ , 

... . „ ,, . . , , r quency, and L„ as the gas surface density (Toomre 1964, 1977 Gol- 

planets could form via direct collapse - fragmentation of a self- n — — r^, £ — „ i .„,_i H= 1 ' . — ■ 

.. .. , , dreich & Lynden-Bell 1965 1 appears to otter some guidance. And 

gravitating region m a proto-planetary disk - as opposed to accre- h — : — r~ ? — ^ „ , ^ .• i r ,n 

. — _. — £J1 . — . — n . indeed, tor Q <C 1 , most of the mass in disks fragments into self- 



tion onto planetesimals (see e.g. Boley et al. 2006; Boley 2009, . '. . '. . , ... 

_ — , „ , , j ,° I — t v __fJI__L gravitating clumps in roughly a single crossing time. But this was 

Dodson-Robinson et al. 2009, Cat et al. 2008, 2010; Vorobyov & ^ . jr , , , j . j, , 

„ . l i ^»,Jl i „ i ^^J. ii ^ 1 1 .u»., i' derived tor a smooth, homogeneous disk, dominated by thermal 

Basu2010 1 Johnsonetal.2010,Boss2011,Stamatellosetal.2011, ... ' J .. . . i .... . 

pressure with no cooling, so does not necessarily imply stability in 



Seager et al. 201 1 , and references therein). Even if this is not a sig- , , ^ n w,.ii j i- • 

; — • — r-r 1 . . . . . . ... , any turbulent system. Gammie (2001 1 studied a more realistic case 

nificant channel tor planet formation, it is clearly critical to under- 



stand the conditions needed to avoid fragmentation. This is espe- 
cially demanding because such fragments need only form an order- 
unity number of times over the millions of dynamical times most 



of a turbulent disk with some idealized cooling, and showed that if 
the cooling time exceeded a couple times the dynamical time, the 
disk could maintain its thermal energy (via dissipation of the turbu- 
lent cascade), with a steady-state Q > 1 and transsonic Mach num- 
bers (powered by local spiral density waves), thus avoiding runway 
catastrophic fragmentation. 

* E-mail:phopkins@caltech.edu However, many subsequent numerical simulations have shown 
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that, although catastrophic, rapid fragmentation may be avoided 
when these criteria are met, simulations with larger volumes (at 
the same resolution) and/or those run for longer timescales still 
eventually form self-gravitating fragments, even with cooling times 
as long as ~ 50 times the dynamical time (see references above 



and e.g. Rice et al.| 2005 


2012 


Meru & Batej|2011b|a| 


Paardekooper et al.|2011| 


3 aarde 


cooper|2012). Increasing 



resolution and better resolution of the turbulent cascade (higher 
Reynolds numbers) appear to exacerbate this - leading to the ques- 
tion of whether there is "true" stability even at infinitely long cool- 
ing times {Meru & Bate 2012). And these simulations are still only 
typically run for a small fraction of the lifetime of such a disk (or in- 
clude only a small fraction of the total disk mass, in shearing-sheet 
models), and can only survey a tiny subsample of the complete set 
of parameter space describing realistic disks. 

1.2 The Role of Turbulent Density Fluctuations 

Clearly, the theory of disk fragmentation requires revision. But al- 
most all analytic theoretical work has assumed that the media of 
interest are homogeneous and steady-state (though see Kratter & 
|Murra y-Clay 201 1 and references above), despite the fact that per- 
haps the most important property of turbulent systems is their in- 
homogeneity. In contrast to the "classical" homogenous models, 
[Paardekooper (2012) went so far as to suggest that fragmentation 
when <2 > 1 may be a fundamentally stochastic process driven by 
individual, random turbulent density fluctuations - and so can never 
"converge" in the formal sense described by the simulations above. 
That said, simulations of idealized turbulent systems in recent years 
have led to important breakthroughs; in particular, the realization 
that compressible, astrophysical turbulence obeys simple inertial- 



range velocity scalings (e.g. 


Ossenkopf & Mac Low|2002| |Feder-| 


|rath et al.|2010||Block et al. 


20101 Bournaud et al.|2010|l, and that 



- at least in isothermal turbulence - the density distribution driven 
by stochastic turbulent fluctuations develops a simple (log-normal) 
shape, with a dispersion that scales in a predictable manner with 
the compressive Mach number ( Vazquez-Semadeni 1994; Padoan 
|et al.|1997||Scalo et al.|1998||Ostriker et al.|1999) . 

Recently, in a series of papers (Hopkins 2 012b| |2011| 
201 2d e c), we showed that the excursion-set formalism could be 
applied to extend these insights from idealized simulations, and 
analytically calculate the statistics of bound objects formed in the 
turbulent density field of the ISM. This is a mathematical formu- 
lation for random-field statistics (i.e. a means to use the power 
spectra of turbulence to predict the statistical real-space structure 
of the density field), well known from cosmological applications as 
a means to calculate halo mass functions and clustering in the "ex- 
tended Press-Schechter" approach from |Bond et al.|fl991| >. This is 
a well-known theoretical tool in the study of large scale structure 
and galaxy formation, and underpins much of our analytic under- 
standing of halo mass functions, clustering, mergers, and accretion 
histories (for a review, see Zentner 2007 1. The application to turbu- 
lent gas therefore represents a means to calculate many interesting 
quantities analytically that normally would require numerical sim- 
ulations. In Hopkins ( 2012b) (hereafter Paper I), we focused on the 
specific question of giant molecular clouds (GMCs) forming in the 
ISM, and considered the simple case of isothermal gas with an ex- 
actly lognormal density distribution. We used this to predict quan- 
tities such as the rate of GMC formation and collapse, their mass 
function, size-mass relations, and correlation functions/clustering, 
and showed that these agreed well with observations. In Hopkins 
( |2012c[ > (hereafter Paper II), we generalized the models to allow ar- 
bitrary turbulent power spectra, different degrees of rotational sup- 



port, non-isothermal gas equations of state, magnetic fields, inter- 
mittent turbulence, and non-Gaussian density distributions; we also 
developed a time-dependent version of the theory, to calculate the 
rate of collapse of self-gravitating "fragments." 

1.3 Paper Overview 

In this paper, we use the theory developed in Paper I & Paper II 
to calculate the statistics of fragmentation events in Keplerian, sub 
and trans-sonically turbulent disks, with a particular focus on the 
question of fragmentation in proto-planetary disks. We develop a 
fully analytic prediction for the probability, per unit mass and time, 
of the formation of self-gravitating fragments (of a given mass) 
in a turbulent disk. In addition to providing critical analytic in- 
sights, this formulation allows us to simultaneously consider an 
enormous dynamic range in spatial, mass, and timescale, and to 
consider extremely rare fluctuations (e.g. fluctuations that might 
occur only once over millions of disk dynamical times), which is 
impossible in current numerical simulations. We will show that this 
predicts "statistical" instability, fragmentation, and the formation 
of candidate "direct collapse" planets/stars even in the "classical" 
stability regime when 2> 1 and cooling is very slow. However, 
these events occur stochastically, separated by much larger aver- 
age timescales than in the catastrophic fragmentation which occurs 
when Q < 1. We will show that this naturally explains the appar- 
ently discrepant simulation results above, and may be important 
over a wide dynamic range of realistic disk properties for forma- 
tion of both rocky and giant planets. 

2 "CLASSICAL" VS. "STATISTICAL" STABILITY 

Motivated by the discussion above, and the models developed in 
this paper, we will introduce a distinction between two types of 
(in)stability: "classical" and "statistical." 

By "classical" instability, we refer to the traditional regime 
where Q < 1, in which small perturbations to the disk generically 
grow rapidly. A "classically unstable" disk will fragment catas- 
trophically, with a large fraction of the gas mass collapsing into 
self-gravitating objects on a few dynamical times. A "classically 
stable" disk will survive many dynamical times in quasi-steady 
state, and gas at the mean density will not be self-gravitating. 

By "statistical" instability, we refer to the regime where an 
inhomogeneous disk experiences sufficiently large fluctuations in 
density such that there is an order-unity probability (integrated over 
the entire volume and lifetime of the disk) of the formation of some 
region which is so over-dense that it can successfully collapse under 
self-gravity. A "statistically stable" disk has a probability much less 
than unity of such an event occurring, even once in its lifetime. 

Classically unstable disks are always statistically unstable, 
and statistically stable disks are always classically stable. However, 
we argue in this paper that there is a large regime of parameter 
space in which disks can be classically stable, but statistically un- 
stable. In this regime, disks can, in principle, evolve for millions 
of dynamical times with g> 1, and nearly all the disk mass will 
be stable, but rare density fluctuations might form an order-unity 
number of isolated self-gravitating "fragments." 

3 MODEL OUTLINE 

Here, we present an order-of-magnitude, qualitative version of the 
calculation which we will develop rigorously below. This serves to 
illustrate some important scalings and physical processes. 

Consider an inhomogeneous disk (around a star of mass M») 
with scale radius r„, scale-height h, mean surface density E and 
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Figure 1. Here we show the rate of formation of self-gravitating clumps as a function of mass. We plot the probability per unit time, per log-interval in mass, 
of the formation of a self-gravitating gas overdensity (turbulent density fluctuation) which can collapse - i.e. a "fragmentation event" and candidate planet 
formation via direct collapse. Mass is in units of p? where M,i is the disk mass and p = /M, is the disk-to-star (central object) mass; the probability is in 
units of fi~~ £7, where fl = 2tt / f 0I bit is the Keplerian frequency. We define a "reference model" (M c = 1, p = 2, Q = \,b = \ /2, 7 = 1, T = 0) and vary each 
parameter in turn. (1) Large-scale compressive (longitudinal) Mach number _M C = bM[h]. Fragmentation is exponentially suppressed when M c <C 1, and 
develops on a single crossing time over a broad mass range when M c 3> 1. There is a characteristic mass ~ p 2 M^. (2) Turbulent spectral index (E(k) oc k~ p ), 
p = 2 is Burgers (highly compressible) turbulence and p = 5/3 is Kolmogorov (incompressible); values outside this range are rare. This has weak effects, 
but shallower spectra give more power on small scales. (3) Global Toomre parameter (2 ~ 1 for marginal classical stability). 2> 1 exponentially suppresses 
(but does not eliminate) fragmentation. 2< 1 leads to "catastrophic" fragmentation in a single crossing time on a broad range of mass scales. (4) Fraction 
b of turbulent velocity in compressive modes (Ai c = bAi). In purely compressive turbulence 6 = 1, pure solenoidal turbulence (and/or cases with strong 
magnetic fields) b = 1 /3, and random forcing b = 1/2. At fixed compressive Aic (not fixed Ai), this has a weak effect. (5) Equation of state polytropic index 
7 (c 2 oc p 1_ '). "Soft" 7 < 1 produce a much broader spectrum of fragmentation on small scales, as compared to isothermal (7 = 1). 7 = 4/3 corresponds to 
a radiation-pressure supported disk with no cooling; 7 = 5/3 to an adiabatic disk with no cooling. These cases suppress fragmentation on smaller scales, but 
still form collapsing regions via turbulent density fluctuations on larger scales with nearly the same integrated probability (normalization). (6) Intermittency 
parameter T (see Paper II; Appendix B). T = is no intermittency; T = 0.05 corresponds to the intermittency model of |She & Lev eque 1 1994), appropriate 
for incompressible Kolmogorov turbulence; T = 0. 12 to the model of Boldyrev ( 2002 1, for more intermittent compressible super-sonic turbulence. 



mid-plane density po ~ E/(2A), and Toomre Q > 1. Although the 
disk is classically stable (and so not self-gravitating at the mean 
density), if a local "patch" of the disk exceeds some sufficiently 
large density then that region will collapse under self-gravity. The 
most unstable wavelengths to self-gravity are of order the disk scale 
height ~ h: roughly speaking, a gas parcel of this size will collapse 
under self-gravity if it exceeds the Roche criterion (overcomes tidal 
forces): p > M,/ri ~ gp [] 

So for Q > 1 we see that gas at the mean density will not col- 
lapse. However, turbulence produces a broad spectrum of density 
fluctuations. For simple, isothermal turbulence with characteristic 
(compressive) Mach number M = (v 2 urb ) 2 /c s , the (volumetric) 
distribution of densities is approximately log-normal, with disper- 
sion a lnp ~ Vln[l+.M 2 ] (« M for M < 1). 

1 To derive this, we assume a Keplerian disk, Q 2 S3 GM„ / r\ , and vertical 
equilibrium, h R3 c s /Q for a disk supported by thermal pressure. 



Since density fluctuations exceeding p > Qpo (ln(p/po) > 
In Q) can collapse, we integrate the tail of the log-normal density 
probability distribution function (PDF) above this critical density to 
estimate the probability P c ~ erfc[(ln(2)/(\/2crinp)], per unit vol- 
ume, of a region being self-gravitating at a given instant. Since we 
are considering regions of size ~ h, the disk contains approximately 
~ {r,/h) 2 independent volumes, and so (assuming P c is small) the 
probability of any one such volume having a coherent volume- 
average density above the critical threshold is P r ~ {r*/h) 2 P c . 

Now consider that a typical proto-planetary disk (at a few AU) 
might have Q ~ 100 and M. ~ 1, with /i/r* ~ 10 _1 . In this case, 
P c ~ 3 x 10~ 8 and P v ~ 3 x 10~ 6 are extremely small! 

However, this analysis applies to the disk viewed at a sin- 
gle instant. Turbulent density fluctuations evolve stochastically in 
time, with a coherence time (on a given scale) about equal to 
the bulk flow crossing time ~ h/v, (since this is the time for 
fluid flows to cross and interact with a new independent region 
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Figure 2. Top: Probability (integrated over all masses in Fig. ^ per 
unit time of a fragmentation event vs. compressive Mach number Aic, 
for other varied parameters. At fixed Ai c and Q, the turbulent power 
spectrum shape (p), forcing mechanisms and presence/absence of mag- 
netic fields (b), equation-of-state (7), and intermittency have small ef- 
fects. In all cases the rate grows very rapidly when M c > 0.1, and is 
suppressed when Q ^ 1. Bottom: The minimum value of Q for statis- 
tical stability, as a function of _M C . For Q > Qmin, the time-integrated 
probability of a fragmentation event Pf rag i nl = J df dA^g/dt is small (here 
we choose P™^ < 0.1, but the exact choice has only a weak logarith- 
mic effect on the result). The three lines span a plausible range for 
the "lifetime" of a turbulent protoplanetary disk (recall, ~yr at 

Jupiter, and typical p, < 0.1). The curves approximately follow Qmin ~ 
0.5 exp (v^lnCTime/^Q- 1 ) \n(\+M 2 )) ~ 0.5 exp (6 ^/\n{\+M 2 )) 
(see Eq. |16|36( . 



of size ~ h). So the density PDF is re-sampled or "refreshed" 
on the timescale ~ h/v, ~ l/(A4Q). But this is just a dynam- 
ical time, which is very short relative to a typical disk lifetime 
(Q- 1 ~yr at ~ 10 AU). If the disk survives for a total timescale 
to, then, the entire volume is re-sampled ~ Mto^l independent 
times. So the probability, integrated over time, of any one of these 
volumes, at any one time, exceeding the self-gravity criterion, is 
< g ~ M to n (V* /hf erfc[ln<2/( V2a lnp )]. 

If a typical disk with parameters above survives for ~ Myr, or 
~ 10 s crossing times, we then obtain a time-and-volume-integrated 
probability Pf"^ g ~ 1 of at least a single stochastic "fragmenta- 
tion event" driven by turbulent density fluctuations. The mass of 
the self-gravitating "fragment" will be ~ (47r/3)/! 3 p ~ 4Qh 3 po 



(~ 0.1 — lMjupim, for these parameters in a minimum mass solar 
nebula). And, despite the fact that the average timescale between 
such events may be long (~Myr), if a fragment forms, it forms 
rapidly (in ~yr) on the turbulent crossing time ~ \/{MQ.) ~yr 
and has a short collapse/free-fall time ~ 1/^/Gp ~ fi _1 ~yr. De- 
spite having Q ~ 100, then, we estimate that this disk could be 
statistically unstable! 

For otherwise fixed h/r, ~ 0.1 and M. ~ 1, P"^ g declines ex- 
ponentially with increasing Q. So for Q 2> 100, <C 1 and such 
a disk is both classically and statistically stable. For Q < 1, the disk 
is classically unstable, and even material at the mean density (i.e. 
an order-unity fraction of the mass) begins collapsing on a single 
dynamical time. But for 1 < Q < 100, such a disk is classically 
stable, but statistically unstable. 

Below, we present a more formal and rigorous derivation of 
statistical (in)stability properties, and consider a range of disk prop- 
erties. But the simple order-of-magnitude arguments above capture 
the most important qualitative behaviors we will discuss. 



4 THE MODEL: TURBULENT DENSITY 
FLUCTUATIONS 

As mentioned above, turbulence in approximately isothermal gas 
(neglecting self-gravity) generically drives the density PDF to a ap- 
proximate lognormal (a normal distribution in In (p)). This follows 
from the central limit theorem (see |Passot & V azquez-Se madeni| 
[1998) |Nordlund & Padoan|[l999) , although there can be some 
corrections due to intermittency and mass conservation ( |Klessen| 
2000, Hopkins 2012a). And in the simplest case of an ideal box 
of driven turbulence, the variance 5 simply scales with the driving- 
scale Mach number M as S = In [1 + b 2 M 2 ] (where b is a constant 
discussed below). These scalings have been confirmed by a huge 
number of numerical experiments with M ~ 0.01 — 100 , sampling 
the PDF down to values as low as ~ 10 -10 in the "tails" {Federrathl 
let al.|2008|[2TJT0l|Federrath & Klessen|20T2l|Schmidt et al.|2009[ 
|Price & Federrath 2010; Konstan dmet al.|2012fr . 

This is true in both sub-sonic and super-sonic turbulence (see 
references above), as well as highly magnetized media (with the 
magnetic fields just modifying b or the effectively compressible 
component of M; s ee|Kowal et al.|2007||Lemaster & Stone|200"9"{ 
|Kritsuk et al.|2011[ |Molina et al.|2012) , and even multi-fluid me- 
dia with de-coupled electrons, ions, neutral species, and dust (see 
|Downes|2012) . And although the density PDF is not exactly log- 
normal when the gas is no longer isothermal, the inviscid Navier- 
Stokes equations show that the local response is invariant under the 
substitution S(M) -> S(M j p) = ln[l +bvj/c 2 (p)], allowing an 
appropriately modified PDF which provides an excellent approxi- 
mation to the results in simulations (see Scalo et al. 1998 ; Pa ssot &| 
Vazquez-Semadeni 1998). 

In Paper I-Paper II, we use these basic results to show how 
excursion set theory can be used to analytically predict the statisti- 
cal structure of turbulent density fluctuations. The details are given 
therein, but we briefly summarize them here. Consider, for sim- 
plicity, the isothermal (lognormal) case: if density fluctuations are 
lognormal, then the variable <5(x) = ln[p(x) / po\ +S/2, where p(x) 
is the density at a point x, po is the global mean density and S is the 
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variance in In p, is normally distributed according to the PDF^] 

1 



Po(S\S) = 



■. exp 



2S 



(2) 



More generally, we can evaluate the field 5(x \ R), which is the 5(x) 
field averaged around the point x with some window function of 
characteristic radius R; this is also normally distributed (see Paper 
II, Appendix F), with a variance at each scale S(R) that is directly 
related to the density power spectrum. 

Paper I derives the "first-crossing" distribution for the gen- 
eral form of these fields. This corresponds to the mass and initial 
size spectrum of regions which are sufficiently dense so as to be 
self-gravitating averaged on the scale R (specifically defined as the 
largest scale on which the region is self-gravitating, i.e. exclud- 
ing bound sub-units already counted "within" the parent, although 
these can be counted separately if desired). This corresponds to the 
statistics of regions where 5(x|i?) > B(R), where B(R) (the "bar- 
rier") is some (scale-dependent) critical value. In Paper I we derive 
S(R) and B(R) from simple theoretical considerations for all scales 
in a galactic disk and/or molecular cloud. However, the derivation 
proceeds almost identically for a proto-planetary disk, following 
the most general form presented in Paper II, which we outline here. 

It is well-established that the contribution to density vari- 
ance from the velocity variance on a given scale goes as S ~ 
ln(l + M 2 ), where M is the Mach number. For a given turbulent 
power spectrum, then, S(R) is determined by summing the contri- 
bution from the velocity variance on all scales R' > R: 



S{R, p) 



\W(k, R)\ 2 ln [ 



b 2 vj(k) 



C]{ P , k)+K?k- 



A\nk (3) 



where W is the window function for the smoothing] Vf (it) is the 
turbulent velocity dispersion averaged on a scale k (trivially related 
to the turbulent power spectrum), c s is the thermal sound speed 
(both c, and 5 can depend locally on p if the gas is not isother- 
mal), and b ~ 1 is the fraction of the turbulent velocity in com- 
pressive (longitudinal) motions (discussed below). Here k is the 
epicyclic frequency; since we are interested in Keplerian disks, we 
take k as Q, the disk orbital frequency. Note that on large scales, 
angular momentum (ft 2 k~ 2 ) enters in a similar way to c 2 and sup- 
presses fluctuations, which follows directly from the form of the 
dispersion relation for density perturbations (e.g. |Lin et al.|1969| 
|Toomre|1977||Lau & B erlin 1978 ); accounting for this is necessary 
to ensure mass conservation. 

Since we are interested in the formation of self-gravitating gas 
objects, we define B(R) corresponding to the critical density aver- 
aged on a given scale, p a u(R), at which an overdensity will col- 
lapse. Given 5(R) = \n[p(R)/p ] + S/2 defined above, then B(R) 
follows from the dispersion relation for a density perturbation in a 
disk with self-gravity, turbulence, thermal and magnetic pressure, 
and angular mo mentum/shear {Va ndervoort 1970; Aok i et al.|1979| 
|Elmegreen| 19871 |Romeo| 1 992): 



B(R) = In 



Pca(R)\ , S(R) 



Pi) 



+ ■ 



(4) 



2 The +S/2 term in S is required so that the integral of pPo(p) correctly 
gives po with (<5) = 0. 

3 For convenience we take this to be a &-space tophat: W = 1 for k < l/R, 
W = otherwise. But we show in Paper I and Paper II (Appendix G) that 
this has little effect on our results. 



where p cr j t is the critical density above which a region is self- 
gravitating. This is the (implicit) solution to 

Peril (g) _ Q L + h 

po 2 ft \ R 



r a 2 {R, p mt ) h , 2 g] 
I a 2 g (h,p ) R ^ hi 



(5) 



where po is the mean midplane density of the disk, h is the 
disk scale height, ft = k/£1 = 1 for a Keplerian disk, and Q = 
(<j g [h, po] ft)/(7rGE gas ) is the Toomre Q parameter. The "total" 
dispersion a g is 



a 2 g (R, p) = c 2 (p) + (v 2 (R)) + vi(p, R) 



(6) 



(va is the Alfven speed). A full derivation is given in Paper I; but we 
stress that this not only implies/requires that a region with p(R) > 
pent (R) is gravitationally self-bound, but also that such regions will 
not be unbound by tidal shear (i.e. have sizes within the Hill radius) 
and that they will not be unbound/destroyed by energy input from 
the turbulent cascade (shocks and viscous heating). 

The map between scale R and the total mass in the collapsing 
region is 



M(R)=4npcnth : 



TR 2 



l + : 



-1 



(7) 



D exp (-?> 

It is easy to see that on small scales, these scalings reduce to the 
Jeans+Hill criteria for a combination of thermal, turbulent, and 
magnetic support, with M = (4n/3) p cn tR 3 ', on large scales it be- 
comes a Toomre-like criterion with M = 7rE cr i t R 2 . 

For any B(R) and S(R), Paper I shows that the instantaneous 
"first-crossing" mass function (i.e. instantaneous mass function of 
collapsing objects, uniquely defined to resolve the "cloud-in-cloud" 
problem) is 



An 
AM 



Peril (M) 

M 



ff{M) 



dS I 
IdM I 



(8) 



where f/(S) is a function shown by Zhang & Hui ( 2006 i to be the 
solution of the Volterra integral equation: 



with 



/ / (5)=g 1 (5)+ / AS'f f (S')g 2 (S,S') 
Jo 

fi(S)=[2| 
g2(S,S')= [- 



dB\ 5(5)' 
~dS\ + ~S~. 
B(S)-B(S' 
S-S' 



Po(B(S)\S) 



dB 
dS 



(9) 



(10) 



(11) 



P [B{S)-B(S')\S'-S] 



In Paper II, we further generalize this to fully time-dependent 
fields. In statistical equilibrium, the density field obeys a modi- 
fied Fokker-Planck equation with different modes evolving stochas- 
tically - they follow a damped random walk with a correla- 
tion time equal to the turbulent crossing time for each spatial 
scale/wavenumber. This allows us to directly calculate the prob- 
ability per unit time of the formation of any bound object in the 
mass function above. The exact solution requires a numerical ap- 
proach described in Paper I (§ 7) and Paper II (§ 9-10), but to very 
good approximation (for Md«/dlnM < po) this is just 



AP(M, At) i 



An At 



AM t{M) 



(12) 



where r(M) = R[M]/v t (R[M]) is the turbulent crossing time on the 
scale R[M] corresponding to the mass (Eq.|7J. This is because the 
coherence time of density fluctuations on a given scale is just the 
crossing time. 



© 0000 RAS, MNRAS 000, 000-000 



6 Hopkins 



Now consider a disk, or disk element (if we consider a series 
of cylindrical annuli at different disk-centric radii) with total mass 
Md- By the definitions used in Paper I-Paper II, this means the total 
"effective volume" is M c \ / po, i.e. that the total (integrated) number 
of objects per unit mass is 



dN 



Md dn 



dlogM po dlogM 



M po I dlogM I 



(13) 



For any polytropic gas, we can factor out all dimensional 
parameters, and work in units of h and po- Mass then has units 
of poh 3 . But since, for a disk in vertical equilibrium h = cr s /f2, 
and for a Keplerian disk n = Q, = (CM,/-, -3 ) 1 ' 2 (where M, is 
the central mass and r» is the disk-centric radius), we have Q = 
(ff,[A]/s)/(jrGE IO ) = hn 2 /(7YGE g , s ) = {h/rjM./faXg.ti). 
So we can write (h/r*) = Qp where 



7rE E . 



M, 



(14) 



where is the disk mass within R. For an exponential vertical 
profile used to define our dispersion relation, E gas = 2 po h, so we 
can use this and (h/r,) — Qp to link poh 3 — (27r) _1 (pQ) 2 M d . 
We can then remove the local quantities po and h and define mass 
in the much simpler global units of p and Mj. The units of dN 
defined above are Mj / (poh 3 ), so this can be similarly re- written. 
And since the disk crossing time scales as ~ h/a g [h] = Q.~\ Q. 
provides a natural time unit. 

Having defined units, the model is completely specified by di- 
mensionless parameters. These are the spectral index p of the tur- 
bulent velocity spectrum, E(k) oc k~ p (v 2 (R) oc R p ~ 1 ), and its nor- 
malization, which we define by the Mach number on large scales 
Ml = (v 2 (h))/c 2 , as well as the Toomre parameter Q. We must 
also specify the parameter b, i.e. the mean fraction of the velocity in 
compressive modes. This is h = 1 for purely compressively forced 
turbulence, b = 1/3 for purely solenoidally forced turbulence, and 
b = 1/2 for random forcing. But we can almost completely fac- 
tor out the dependence on this parameter if we simply define the 
compressive component of the turbulence, i.e. work in units of the 
compressive Mach number M c = bM. 

In Paper II, we show how these equations generalize for the 
cases with non-isothermal gas, intermittent turbulence, and non- 
isotropic magnetic fields. The qualitative scalings are similar, but 
the math becomes considerably more complicated (and the "first- 
crossing distribution" //(M) and its evolution in time must be 
solved via a numerical Monte Carlo method rather than analyti- 
cally). The presence of intermittency makes the density PDF non- 
Gaussian in the "tails," but this can be entirely encapsulated in a 
modified form of Eq. [9] above (see Paper II, Appendix D), leaving 
the rest of our derivation intact. For non-isothermal cases, we must 
replace c s — > c s (p) (in calculating both the variance and critical 
density), and again allow for the density PDF to be non-Gaussian 
(with a skew towards lower/higher densities as the equation of state 
is made more or less "stiff," respectively; see § 3-4 in Paper II). 
Modulo these changes, however, our derivations (and the changes 
made to specify to the Keplerian disk case) are identical for the 
simple case of a polytropic gas where c 2 oc p 7_1 (7 is the poly- 
tropic index). Magnetic fields can produce global anisotropy, but 
this can be simply absorbed into the form of //(M) as well; their 
largest effect in suppressing fluctuations actually manifests as a 
field-strength-dependent value b (Paper II, § 6, and e.g. |Kowal et al.| 
|2007l |Lemaster & Stone|2009[|Molina et al.|2012| >. For the strong- 
field limit, however, this is just similar to the pure-solenoidal turbu- 
lence case (in both cases, there is only one spatial dimension along 



which compression is possible); so is within the range of b varia- 
tions we will consider. 



5 RESULTS: FRAGMENTATION RATES AND MASS 
SPECTRA IN THE GENERAL CASE 

In Fig. [T] we now use this to predict the mass spectrum - specifi- 
cally the probability per unit time, per unit mass - of the formation 
of self-gravitating regions, as a function of various properties of the 
systemfl 

We define our "reference" model to be isothermal (7 = 1), 
non-intermittent (i.e. log-normal density PDF), with turbulent spec- 
tral slope p — 2 (appropriate for compressible turbulence), Toomre 
2 = 1, and have b = 1/2 (random forcing) with rms three- 
dimensional compressive Mach number on the largest scales M. c = 
1. But we then vary all of these parameters. For now, we treat each 
as free - in other words, we make no assumptions about the specific 
mass profile, temperature structure, cooling chemistry, or other mi- 
crophysics of the disk. These microphysics are, of course, what ulti- 
mately determine the value of the model parameters (and may build 
in some intrinsic correlations between them); but for any given set 
of parameters in Fig. [T| the prediction is independent of how the 
microphysics produce those parameters. 

In general, the shape of the mass spectrum is similar re- 
gardless of these variations. This is peaked (very approximately 
lognormal-like) around a characteristic mass ~ 0.1 - \p 1 M i . For 
a relatively massive disk Md ~ 0.1M*, then, we expect character- 
istic masses in such events of ~ 10~ 3 M*, corresponding to gas 
giants. If however the disk mass is lower, Mj ~ 0.0 1M„, then this 
becomes ~ 10 -6 M„, typical of rocky (Earth-mass) planets. 

The dependence on the slope of the turbulent spectrum is quite 
weak: Kolmogorov-like (p = 5 /3, appropriate for incompressible 
turbulence) spectra give nearly identical results; shallower spectra 
slightly broaden the mass range predicted (since M c declines more 
slowly at small scales), but these are not seen in realistic astrophys- 
ical contexts. Likewise, at fixed A4 C , there is some dependence on 
b (because of how M — M c /b enters into the critical density for 
collapse), but this is small in this regime, because turbulence is not 
the dominant source of "support" resisting collapse. The effects of 
intermittency are very weak, since they only subtly modify the den- 
sity PDF (see Paper II). Even changing the equation of state has 
surprisingly mild effects. A "stiffer" (higher-7) equation of state is 
more resistive to fragmentation on small scales and leads to a more 
sharply peaked spectrum. But for fixed A4 C , the variance near the 
"core" of the density PDF is similar independent of 7, so this does 
not have much effect on the normalization of the probability distri- 
bution. We stress that the medium could have a perfectly adiabatic 
equation of state with no cooling, and ;/ it had the plotted Mach 
number and Toomre Q, our result would be identical (and the prob- 
ability of fragmentation would still be finite). Very soft equations 
of state, on the other hand, can lead to a runaway tail of small- 
scale fragmentation, but this is not likely to be relevant. Clearly, 
the largest effects come from varying M c and Q; at A4 C > 1 or 
Q < 1 the mass distribution rapidly becomes more broad, while at 



4 For simplicity, we will refer to these calculations as if the disk is constant 
surface-density out to some maximum radius, with total mass Mj and Q (the 
orbital velocity) defined at that radius, and constant Toomre Q parameter. 
However, the results shown could be applied to any radius in a disk with 
any mass profile (provided it is still Keplerian), with Mj defined as the disk 
mass inside a radial annulus, and Q and Q evaluated at that same radius. 
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M c ^ 1 or Q > 1 the characteristic mass remains fixed but the nor- 
malization of the probability becomes exponentially suppressed. 

This is summarized in Fig.[2] Since the mass range of expected 
"events" is relatively narrow, we integrate over mass to obtain the 
total probability of an event per unit time, 



df 



dlogM 



AN 



dlogMdf 



(15) 



and plot this as a function of M c for different parameter choices. 
As before, most parameters make a surprisingly small difference; 
M c and Q dominate. 

5.1 A General Statistical Stability Criterion 

Formally d/Vf rag /df is always non-zero for M c > 0; but we see that 
for for M c < 1/2, the probability per unit time of forming a self- 
gravitating fluctuation drops rapidly. However, recall that the total 
lifetime of e.g. a proto-planetary disk is many, many disk dynami- 
cal times ~ Q, ~ 1 , and our "time unit" is fi 2 Q ~ 1 . Consider a typical 
lifetime of TMyr = Tdisk/Myr ~ 1; then the disk at ~ 10 AU around 
a solar-mass star (where fl « lyr -1 ) experiences ~ 10 6 dynamical 
times. For a disk-to-total mass ratio of fi ~ 0.1, this is ~ 10 8 "time 
units," so if we integrate the probability of a fragmentation event 
over the lifetime of the disk we obtain an order-unity probability 
even for d/Vf rag /dr ~ 10~ 8 in the units here (i.e. M c as small as 
~ 0.15). Of course, following such an integration in detail requires 
knowing the evolution of the disk mass, Q, M c , etc. But we can ob- 
tain some estimate of the value of Q required for statistical stabil- 
ity (ensuring the probability of fragmentation events is negligible) 
by simply assuming all quantities are constant and integrating over 
an approximate timescale, also shown in Fig. [2] Here we take our 
standard model, and consider three timescales in units of /i 2 Q~ 1 (a 
factor of ~ 100 shorter and longer than the value motivated above), 
and consider the minimum Q at each M c needed to ensure that the 
time-integrated probability of a fragmentation event is <C 1 . 

This minimum Q mm is ~ 1 at M ~ 0.1; equivalently, disks 
with 2 ~ 1 and M c > 0. 1 have an order-unity probability of at least 
one stochastic fragmentation event over their lifetime (P^L ~ 1). At 
larger _M r > 0.3 - 0.5, g > 3 - 5 is required for fg£ g < 1 ; and by 
sonic Mach numbers M c ~ 1—3, <2 > 40 — 1000 is required for 
< g «l- 

We can approximate the scaling of Q m m(M c ) by the follow- 
ing: recall that the critical density near the Toomre scale h scales 
approximately as In (pmt/po) ~ In (2 Q) (Eq.l^l, while the density 
dispersion scales as a\ np ~ \J (1 + M 2 ) (Eq.pL As noted above, 
if the system evolves for a total timescale to = fo/(//f2 _1 ) (time 
in our dimensionless units), then an event with probability per unit 
time P ~ 1 /to has an order-unity probability of occurring. If the 
probabilities are approximately normally distributed then this is 
just exp(— B 2 / (25)) ~ 1 /to, where B is the barrier and S the vari- 
ance. Since the mass function is peaked near the Toomre scale we 
can approximate both by their values near the "driving scale" ~ h, 
B ss In (2 Q) and 5 w ln(l +Mf). Thus, statistical stability over 
some timescale of interest to requires a Q mm in Fig.|2]of 

i 0.5 exp [ ^2 In {to p~ 2 In ( 1 + M 2 )] 



Qn 



(16) 



For typical values of to, Q, and p - i.e. the typical number of 
independent realizations (in both time and space) of the turbulent 
field in a protoplanetary disk, this becomes 



gmm ~ 0.5 exp (6 y/\n{l+M 2 )) 



(17) 



In other words, a ~ 5 — 6ai np event has order-unity chance of oc- 
curing once over the disk lifetime, so for any M c this implies a 
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Figure 3. Characteristic mass at collapse - i.e. "seed" or "collapse" mass - 
as a function of radius, for a protoplanetary disk with surface density profile 
S = So looo lOOOgcm - 2 (r*/au) -Q , around a star with mass M„ and disk 
temperature calculated including illumination and accretion as described in 
§UJ Specifically, the mass is the mean (M) of the predicted MF as in Fig.JT] 
calculated with Q for this temperature and S, and turbulent b = 1 /2, T = 
0, 7 = 7/5, p = 2, Mc = 0.1 (these parameters have little effect on the 
prediction). The MMSN corresponds to So. 1000 = 1, a = 1 -5, M t = Mq . 
Units are Jupiter masses. The characteristic mass scales as ~ pr Q 2 Mj(< 
r„). This (on average) increases with r t , spanning an Earth-to- Jupiter mass 
range (M Earth = 0.003Mj). 



minimum Q needed to ensure statistical stability in such an extreme 
event. 



6 HOW IS THE TURBULENCE POWERED? 

STATISTICAL STABILITY IN SPECIFIC MODELS 
FOR TURBULENCE 

Thus far, we have considered the general case, varying the Mach 
numbers A4 C independent of other disk properties such as Q, E, 
and 7. However, in a realistic physical model, the mechanisms that 
drive turbulence may be specifically tied to these properties. More- 
over, there may be certain characteristic Mach numbers expected or 
ruled out. In this section, we therefore consider some well-studied 
physical scenarios for the driving of turbulence in Keplerian disks, 
and examine their implications for the "statistical stability" we have 
described above. 

6.1 The "Gravito-turbulent" Regime (Gravity-Driven 
Turbulence) 

Much of the work studying fragmentation in Keplerian disks has 
considered disks with a (locally) constant-cooling rate ("£" disks, 
with t coo i = £f2 - ' locally fixed). In particular, this includes the sce- 
nario of a "gravito-turbulent" steady-state from Ga mmie| ( |2001| >, 
with local instabilities (spiral waves) powering turbulence which 
contributes an effective viscosity and maintains a steady tempera- 
ture and Q ~ 1 . The theory we present above is more general than 
this: we make no assumption about the detailed cooling physics, 
or that the disk is an a-disk, and allow the various parameters Q, 
Mc, 7, etc. to freely vary, whereas many of these are explicitly 
linked in these models. But in the theory above we cannot predict 
these quantities (or their co-dependencies); therefore, this model 
provides a simple and useful way to relate and predict some of the 
otherwise independently free parameters of the more general case, 
and is worth considering in detail. 
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6.1.1 General Settlings 

Consider a cooling rate which is uniform over a region (annulus) of 
the disk, 



^cool — C 



(18) 



If dissipation of gravitational instabilities (e.g. spiral waves) pro- 
vides a source of heating balancing cooling, and £ > 3, the system 
can develop a quasi steady-state angular momentum transport and 
Toomre Q parameter, as in a Shakura & Sunyaev] {1973} a-disk. 
Pringle ( 1981 1 and Gammie ( 2001} showed that in this equilibrium, 
the "effective" viscosity parameter a 



(19) 



9 7 (7-l)C 



is approximately constant. This corresponds to the amplitude of 
density waves 5E/E oc a 1 ^ 2 , leading to a maximum a max « 0.06 
above (hence minimum £ m i„ below) which the local Q < 1 and 
catastrophic fragmentation will occur (see |Gam mie 200 lj |Rice| 
|etll]|2005l |Cossins et af1|2009l >. In an a-disk, the inflow rate is 
also determined, as 



M = 3irac, E gas f2 



(20) 



so this also corresponds to a maximum "classically stable" inflow 
rate below which catastrophic fragmentation will not occur. 

Implicitly, the relations above also define a steady-state Mach 
number. Recall, the dissipation of spiral instabilities is ultimately 
governed by the turbulent cascade. Since the turbulent dissipa- 
tion rate is constant over scale in a Kolmogorov cascade, we can 
take it at the top level, dE/dA df = (p — 1 ) ~~ 1 E gas v 2 f2 (where here 
we consider the rate per unit area) and equate it to the cooling 
rate = [7(7- 1)] _1 EgascJ t^, giving M 2 C ~ (3/2) a. Equiva- 
lently we could have equated the Reynolds stress that leads to a, 
a = (dlnn/dlnr,) -1 7R e y/(E gas c 2 ) with r Rcy = (E gas Sv r 5v^) ; we 
obtain 



1 



(21) 



2 3 7 ( 7 -l)C 

again^We now see how this relates to the theory developed in this 
paper. 

In the language here, increasing £ enters our theory by - as 
discussed in the text - changing the equilibrium balance between 
thermal and turbulent energy, i.e. the Mach number. As cooling be- 
comes less efficient, maintaining the same Q ~ 1 needed to power 
the turbulence (since spiral waves are still generated) requires a 
smaller turbulent dispersion, hence makes the system "more sta- 
ble." 

6.7.2 A Statistical Stability Criterion 

But this also suggests an improved statistical stability criterion, ac- 
counting not just for regions where Q < 1 but also for stochastic 
local turbulent density fluctuations. For a given Q, we have calcu- 
lated (Fig. [2} the Mach number M c above which the system will 
be probabilistically likely to fragment on a given timescale. Eq.|21| 
allows us to translate this to a minimum (. We can do this exactly 
by simply reading off the numerically calculated values, but we 



One subtlety here is that the hydrodynamic Reynolds stress, and/or the 
dissipation on small scales, is dominated by the longitudinal (compressive) 
component. So this relation actually determines _M C , not necessarily M. 
But for our purposes, this is particularly convenient, as it allows us to drop 
the b term from our earlier derivation. 



can also obtain an accurate analytic approximation by the follow- 
ing (see § [8}. Recall that for a system which evolves for a total 
timescale to = to/ (p, 2 (time in the dimensionless units we 
adopt), we obtain the approximate Eq. [T6] for the Q m \ n needed to 
ensure 7^4 « 1: 

Q > \ exp [ s/l ln(ro)ln(l+A12)] (22) 

We can invert this to find the maximum M c for statistical stability 
at a given Q: 



,: .._fDn(2fi)] 2 



M c < exp 



[ Lin (2^ 
L 21m 



1 : 



[ln(2g)] 2 



(23) 



ito J 21nro 
Where the second equality follows from the fact that (for the sys- 
tems of interest) 2 In (to) 2> In (2 Q) is almost always true. 
Combining this with Eq.|2T] we obtain: 

4 In (r ) 



(24) 



(25) 



37(7-1) [ln(2g)p 

__4 ln[f /(/i 2 £r')] 

37(7-1) [ln(2G)] 2 

We can immediately see some important consequences. Be- 
cause of the stochastic nature of turbulent density fluctuations, Cmin 
will never converge in time integration (assuming the disk can 
maintain steady-state mean parameters) - there is always a finite 
(although possibly extremely small) probability of a strong shock 
or convergent flow forming a region which will collapse rapidly. 
However, the divergence in time is slow (logarithmic). The criti- 
cal £ also scales with 7 just as in the Gammie (2001) case; as the 
equation of state is made "stiffer," the Mach numbers and density 
fluctuations are suppressed so faster cooling (lower Q can be al- 
lowed without fragmentation. And £ scales inversely with log (Q), 
so indeed higher-g disks are "more stable," but there is no "hard" 
cutoff at a specific Q value. 

We can turn this around, and estimate the typical timescale for 
the formation of an order-unity number of fragments at a given (, 
obtaining 

(t(N fme ~ 1)) « //JT 1 exp C7(7- 1) Dn(2(2)] 2 ) (26) 

As expected, this quickly becomes large for modest £ and/or Q. 
We stress, however, that this is a probabilistic statement. Although 
the mean timescale between fragment formation events might be 
millions of dynamical times, if and when individual fragments form 
meeting the criteria in the text, they do so rapidly - on of order a 
single crossing time. 

Given our derivation of £ m i n , what do we expect in realistic 
systems such as proto-planetary disks? For a physical disk with 
( > 1 we should expect 7 ~ 7/5 — 5/3 (depending on the gas 
phase), and in equilibrium Q ~ 1; and we should integrate over 
the entire lifetime of the disk, to = r tot ~ 10 6 - 10 10 (with values 
motivated in §|5}. We then expect 

51 l+0.051n(T nt /10 8 ) 



min ' 



(27) 



7(7-1) (l + 1.41n<2) 2 

This is fairly sensitive to Q - note ~ 13/ [7 (7 - 1)] if Q = 
2 instead - and weakly sensitive to to for reasonable variations. 
But this implies that only disks with extremely slow cooling, t coo \ > 
50 fi -1 , (corresponding to steady-state Mach numbers M c < 0.02) 
are truly statistically stable with Q ~ 1 over such a long lifetime. 

6.1.3 Comparison with Simulations 

Now consider the parameter choices in some examples that have 
been simulated. Gammie (2001) considered the case with 7 = 2, 



© 0000 RAS, MNRAS 000, 000-000 



Gravo -Turbulent Planet Formation 9 



and a steady-state Q m 2.46, evolving their simulations for typ- 
ical timescales to ~ 50 SI - ' (though they consider some longer- 
scale runs). Because these were two-dimensional shearing-sheet 
simulations, the appropriate p is somewhat ambiguous, but recall 
{h/r*) = Qp by our definitions, and for the assumptions in Gam- 
( 2001 1 their "standard" simulation corresponds to an h/r* ss 



0.01 Q (where we equate the "full disk size" to the area of the box 
simulated). Plugging in these values, then, we predict £ m i„ = 3.4, in 
excellent agreement with the value £ ss 3 found by trial of several 
values therein. Paardekooper (2012) considered very similar sim- 
ulations but with Q » 1 and all sheets run for t ~ 1000 fT 1 ; for 
this system we predict ( mm = 20.4 - again almost exactly their es- 
timated "fragmentation boundary." |Meru & Bate (2012) and Rice 
et al. (2012) consider three-dimensional global simulations; here 
p = Md/M* =0.1 is well-defined, a more realistic 7 = 5/3 is 
adopted, and the disks self-regulate at Q w 1 ; the simulations are 
run for a shorter time ~ 50— 100 (SI) -1 (where for convenience 
we defined (fi) at the effective radius of the disk, since it is radius- 
dependent, but this is where the mass is concentrated), giving 
Cmin ~ 21 — 25, in very good agreement with where both simula- 
tions appear to converge (using either SPH or grid-based methods). 
This is also in good agreement with the earlier simulations in |Rice| 



et al. 1 2005 1, for 7 = 5/3 (predicting £ m i n = 7.5, vs. their estimated 
6 — 7) and 7 = 7/5 (predicting £ m j„ = 14, vs. their estimated 13). 
Of course, we should naturally expect some variation with respect 
to the predictions, since this is a stochastic process, but we do not 
find any highly discrepant results. 

We should also note that convergence in the total fragmen- 
tation rate in simulations - over any timescale - requires resolv- 
ing the full fragmentation mass distribution in Fig.[T] Unlike time- 
resolution above, this is possible because there is clearly a lower 
"cutoff" in the mass functions (they are not divergent to small 
mass), but requires a mass resolution of ~ 0.01 — 0.1 pf ' Md (de- 
pending on the exact parameters). This is equivalent to a spa- 
tial resolution of e r ~ 0.02 - 0.2Q [/2 h, i.e. a small fraction of 
the disk scale-height h. This also agrees quite well with the spa- 
tial/numerical resolution where (at fixed time evolution) many of 
the studies above begin to see some convergence (e.g. |Meru & Bate| 
|2012|[Rice et al.|2012} , but it is an extremely demanding criterion. 

6.2 The Magneto-Rotational Regime 

In the regime where the disk is magnetized and ionized, the 
magneto-rotational instability (MRI) can develop, driving turbu- 
lence even if the cooling rate is low and Q 2> 1 . We therefore next 
consider the simple case where there is no gravo-turbulent instabil- 
ity fcooi — > oo), but the MRI is present. 

6.2.1 General Settlings 

Given MRI and no other driver of turbulence, Alfven waves will 
drive turbulence in the gas to a similar power spectrum to the hy- 
drodynamic case (within the range we examine where the power 
spectrum shape makes little difference), with driving-scale rms 
{v}) 1 / 2 f& va. In terms of the traditional /3 parameter (ratio of ther- 
mal pressure to magnetic energy density; /3 — > as magnetic field 
strengths increase), /3 = 2c 2 s /v\, so the rms driving-scale Mach 
number is Ai ~ \/2f3~ l . Magnetically-driven turbulence is close 
to purely solenoidal, so b « 1/3 andX c = bM ~ ( V2/3) /3~' /2 . 

We stress, though, that what is important is the saturated lo- 
cal plasma (3 = /3 sat , which can be very different from the ini- 
tial mean field /3o threading the disk. As the MRI develops, the 
plasma field strength increases until it saturates in the fully non- 
linear mode. Direct simulations have shown that for initial fields 
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A> < 10 4 , /3 sat ~ 1/3-2/3 (see e.g. |Bai & Stone|2012[[Fromang 
|et al.|2012] and references therein). The saturation occurs in rough 
equipartition with the thermal and kinetic energy densities - i.e. the 
turbulence is trans-sonic or even super-sonic (M c ~ \/2/3 ~ 0.8). 
In the weak-field limit, however, with /3o > 10 4 , the saturation is 
much weaker, with /3 silt ~ 10 — 20, so M c ~ 0. 1 . 

These same simulations allow us to directly check our simple 
scaling with vaJ the authors directly measure the rms standard de- 
viation in the (linear) density 8 = p/{p), which for a lognormal 
density distribution is (by our definitions) identical to A4 C - For 
four simulations with /3q = 10 2 , 10 3 , 10 3 (but higher-resolution), 
10 4 , they see (midplane) saturation /3 sat = 0.4, 1.1,0.7, 18 and 
(p 2 ) 1/2 /{p) = 0.60,0.43,0.53,0.13 (compared to a predicted 
(p 2 ) [/2 /{p) =M C = (V2/3)P~' /2 = 0.72,0.44,0.55,0.11, re- 
spectively). Moreover, these and a number of additional simula- 
tions have explicitly confirmed that our lognormal assumption (in 
the isothermal case) remains a good approximation for the shape of 
the density PDF (s ee|Kowal et al.|2007[|Lemaster & Stone|2009[ 
|Kritsuk et al.|2011|[Molina et al.|2012) . So for a given M c and Q, 
our previous derivations remain valid. 

6.2.2 A Statistical Stability Criterion 

The strong-field limit therefore leads to large density fluctuations. 
However, strong magnetic fields will also provide support against 
gravity, modifying the collapse criterion; this appears in Eq.[6] But 
for a given /3, this simply amounts (to lowest order) to the re- 
placement c 2 — s> c 2 + v\ = c 2 (1 + 2/3 -1 ). Because Q oc c s , near 
the Toomre scale, this is approximately equivalent to raising the 
stability parameter as Q — > £?cff = <2(v A = 0) +2/3 -1 (where 
Q(va = 0) is the Q including only thermal support). The energy 
and momentum of the bulk flows in the gas turbulence also provides 
support against collapse, so the "effective dispersion" in Eq. [^in- 
cludes all three effects; however this is already explicitly accounted 
for in our previous calculations for any A4. But while the effective 
Q e ff increases in the strong-field limit with /3~ 1,/2 , so does M c , and 
the Q needed for statistical stability on long timescales (Fig. [2J in- 
creases exponentially with M c - so the net effect of MRI is always 
to increase the probability of stochastic collapse. 

Putting this into our general criterion Eq.[l6] we can write the 
statistical stability requirement 



2(va = 0) ; 
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where to = to p~ 2 S7 as before and the latter uses the fact that ft 
is not extremely small in the cases of interest. Integrated over the 
lifetime of the disk, this becomes 
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(29) 



This increases rapidly with increasing magnetic field strength: 
Q mto (vA = 0)wl5 ) 7,3, 1.2 for /3 sal = 1/3,2/3,2,10. 

So MRI with saturation /3 sat < 10 ("seed" /3 < 10 4 ) will make 
even Q > 1 disks statistically unstable, without the need for any 
other source of turbulence. On the other hand, weak-field MRI with 
/3 sa t ^ 10 produces only small corrections to statistical stability. 
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6.3 Convective Disks 

A number of calculations have also shown that proto-planetary 
disks are convectively unstable over a range of radii ( |Boss|2004| 
|Boley et aT1|2006[ |Mayer et aL||2007| >. Most simulations which 
see convection have also seen fragmentation, which has been in- 
terpreted as a consequence of convection enhancing the cooling 
rates until they satisfy thejGammie (2001 1 criterion for fragmenta- 
tion. But |Rafikov|poi)7) and others jCai et al.|200l^|2008] > have ar- 
gued that while convection can and should develop in these circum- 
stances, the radiative timescales at the photosphere push the cooling 
time above this threshold. However, as we have discussed above, 
that would not rule out rarer, stochastic direct collapse events. 

Consider a polytropic thin disk; this is convectively unstable 
when it satisfies the Schwarzschild criterion 



dlnT 7 — 
dlnP > ~Y 



(30) 



Followin g |Lin & PapaloTzoul (1980) , |Bell & Lin| ( [1994) , and 
Rafikovl p0077 " using the fact that disk opacities can be approxi- 
mated by k ~ K()P" r a (P the midplane pressure), this can also be 
writted (l+a)/(4 — 0) > (7 — l)/7- For the appropriate physi- 
cal values, this implies strong convective instability in disks with 
T < 150K (where k is dominated by ice grains) and at higher tem- 
peratures > 1.5 x 10 3 K when grains sublimate, and marginal con- 
vective instability in between. So this should be a common process. 

A convective disk can then accelerate gas via buoyancy at a 
rate comparable to the gravitational acceleration, implying mach 
numbers M 1 w 0.25 — 1 (depending on the driving gradients; recall 
also this is the three-dimensional M) at the scale height where the 
disk becomes optically tiling] Buoyancy-driven turbulence is pri- 
marily solenoidal forcing, so b « 1 /3 while A4 ~ 0.5 — 1, leading 
to a "maximal" A4 C ~ 0.2 — 0.3 (assuming the convection cannot 
become supersonic; this is approximately what is measured in these 
simulations). If this saturation level is independent of Q (provided 
the disk is convectively unstable at all), we then simply need to ex- 
amine Fig.|2]to determine Qmm for statistical stability; from Eq.|16| 
this is approximately 



Q > 0.5 exp [M c \/21nro] ~ 3 



(31) 



Thus, while this does not dramatically alter the behavior of the sta- 
bility criterion Q, it does systematically increase the threshold Q 
for statistical stability by a non-trivial factor. And indeed, in the 
simulations of Mayer et al. ( 2007 1 and Boss ( 2004 1, fragmentation 
occurs when convection is present at radii where <2« 1.4 — 1.8 > 1. 

6.4 Additional Sources of Turbulence 

There are many additional processes that may drive turbulence in 
proto-planetary and other Keplerian disks, but under most regimes 
they are less significant for our calculation here. 

In the midplane of a protoplanetary disk, where large grains 
and boulders settle and are only weakly aerodynamically coupled 
to the gas, Kelvin-Helmholtz and streaming instabilities generate 



6 From mixing-length theory, we can equate the convective energy flux at 
the scale-height F com = 2pC p 7*v 3 /Tiggra-v (where at ~ h, the acceleration 
ggrav ~ £l~h, h m Cj/£7, and for the relevant parameters C p Ri 1.25 X 10 8 
in cgs) to the cooling flux <r SB 7i. This gives us the approximate estimate 
M w 0.5 (r eff /200K) (Sgas/lOOOgcm- 2 )- 1 ^ (fJ-'/yr) 1 / 3 . This agrees 
well with the simulations in the text when M < 1 , but extrapolates to super- 
sonic values at low S and/or large r, , so convection could be considerably 
more important than we estimate if it does not saturate at velocities ~ c s . 



turbulence. However these only operate in a thin dust layer, and ap- 
pear to drive rather small Mach numbers in the gas, so are unlikely 
to be relevant for direct collapse in the gas and we do not consider 
them further (see e.g. |Bai & Stone|2010||Shi & Chiang|2012) . They 
may, however, be critical for self-gravity of those grains participat- 
ing in the instabilities themselves - a more detailed investigation of 
this possibility is outside the scope of this work (since our deriva- 
tion does not apply to a weakly coupled, nearly-collisionless grain 
population), but extremely interesting for future study. 

Radiative instabilities should also operate if the disk is sup- 
ported by radiation pressure, and/or in the surface layer if a wind is 
being radiatively accelerated off the disk by central illumination. 
The former case is not expected in the physical disk parameter 
space we consider; but if it were so, convective, magneto-rotational, 
and photon-bubble instabilities are also likely to be present i |Blaes| 

6 Socrates 2001; Thompson 2008), which will drive turbulence 
that saturates in equipartition between magnetic, radiation, and 
turbulent energy densities, i.e. produce the equivalent of A4 ~ 1 
throughout the disk (giving results broadly similar to the strong- 
field MRI case). In the wind case, the Mach numbers involved can 
be quite large (since material is accelerated to the escape velocity), 
but unless the surface layer includes a large fraction of the mass, it 
is unlikely to be important to the process of direct collapse. 

In the case of an AGN accretion disk, local feedback from 
stars in the disk may also drive turbulence (as it does in galactic 
disks), and this can certainly be significant in the outer regions of 
the disk where star formation occurs (see Thom pson et al.|2005| >. 
In that case the turbulence may even be super-sonic, in which case 
a more appropriate model is that developed in Paper I-Paper II. In 
the inner parts of the disk, though, where the turbulence is sub- 
sonic, we are not necessarily interested in rare single star-formation 
events. 

7 EXAMPLE: PROTOPLANETARY DISK AND 
FRAGMENTATION RADII 

We now apply the statistical stability criteria derived above to a 
specific model of a proto-planetary disk. This is highly simplified, 
but it allows us to estimate physically reasonable sound speeds, 
cooling times, and other parameters, so allows us to ask whether 
our revised statistical stability criteria are, in practice, important. 

7.1 Disk Model Parameters 

For convenience, consider a disk with a simple power-law surface 
density profile 



E = Eo, 1000 1000 gem 2 (r*/au) 



(32) 



The minimum mass solar nebula (MMSN) corresponds to Eo, 1000 ~ 
1 and q ~ 1.5, but we consider a range in these parameters below. 

For a passive flared disk irradiated by a central star with radius 
= and temperature = T„, the effective temperature is (Chiang & 

|Goldreich|19"97l >: 



Teff. • 



1/4 



(33) 



where for a solar-type star a w 0.005 (r*/au) ' + 0.05 (r*/au) 2//7 , 
R, — Rq and T, — 6000 A*. If, instead of irradiation, disk heating 
is dominated by energy from steady-state accretion with some M, 
energy balance requires an effective temperature (i.e. disk flux) 



7cff . 1 



r 3 Mtfi 1 ' 4 

L87T OSB J 



(34) 
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P(frag) « 1 
] P(frag) ~ 1 



1000 



* 100 



lOOOr 



£ 100 




MMSN 



P(frag) « 1 r — r\ P(frag) ~ 1 



Figure 4. Shaded regions show the range of disk temperatures, at a given 
radius around a solar-type star, in which a proto-planetary disk is statis- 
tically unstable (i.e. has an order-unity probability P£j„ ~ 1 of at least 
one "fragmentation" or direct collapse event in a timescale 1 Myr). Top: 
Disk with E = lOOOgcm -3 (r, /au) -1 . Bottom: 10* higher E (E . 1000 = 
10 4 gcm~ 2 ; M,ji s i;(< lOau) ~ 0.05M,). Black is the standard Toomre 
Q < 1 (catastrophic fragmentation). Other shaded regions correspond to 
different mechanisms driving turbulence, with A4 C , Q, etc. calculated self- 
consistently for S(r»), T, and M* (see §[6j. Green: Temperature where 
^fra ~ 1 if the disk is convectively unstable (§ 6.3 Eq. 31), Solid/dashed 
lines correspond to the higher/lower Ai estimatecTTrom convective driving 
in simulations. Pink: Temperature where Pf" a ' g ~ 1 if the disk has a saturated 
(strong-field) magneto-rotational instability (MRI; § |6.2| Eqn |28) ; again 
solid/dashed correspond to stronger/weaker limits on saturation (ft a t = 
0.45 - 0.85, from seed /3 = 10 2 - 10 3 ). Blue: Gravito-turbulence (§[6T| 
Eqn. |24) ; here, higher T corresponds to faster cooling, hence higher M c 
(Eqn|2l|, and increased P^L. Red lines show the calculated T for a disk 
with the given E(r«) from illumination by a solar-type star (dot-dashed) & 
illumination plus accretion with M = 3x 10 — 7 Mq yr — 1 (dotted). Even in 
MMSN, radii > a few au are statistically unstable via gravito-turbulence; 
smaller radii are statistically unstable for M > 3 X 1 0~ 7 Mq yr _ 1 . Strong- 
field MRI is also capable of generating sufficient fluctuations for direct col- 
lapse down to ~ 0. 1 — 1 au. 



where <tsb is the Stefan-Boltzmann constant. The temperature of in- 
terest for our purposes, however, is the midplane temperature r m id, 
since this is where the disk densities are largest and what provides 
the c s resisting collapse; this is related to T^s by a function of opac- 
ity and E which we detail in Appendix [A] But having determined 
T c ff and E, it is straightforward to calculate T m a; the sound speed is 
c s = \J kg T m id I fi where kg is the Boltzmann constant and fj, is the 
mean molecular weight. 




Figure 5. Shaded regions show the range of surface densities which are 
statistically unstable (have an order-unity probability F^„ ~ 1 of a di- 
rect collapse event in a timescale 1 Myr, as Fig.|4j, in a proto-planetary 
disk at a given radius around a solar-type star. For each E(r, ), we cal- 
culate T self-consistently including illumination and accretion with M = 
3x lO _7 M0yr — 1 ((op) or illumination only (bottom); see §|vj Each shaded 
range corresponds to a different candidate source of turbulent density fluctu- 
ations, as in Fig. [4] Red line shows the MMSN, for comparison. Strong MRI 
can promote collapse in a MMSN at all radii; gravito-turbulence in even 
lower-density disks at > au, and smaller radii if M > 3 X 1O _7 M0 yr — 1 . 
If these are not active, convection can promote collapse in higher-density 
disks with E > IOEmmsn- 



This is sufficient to specify most of the parameters of interest. 
In the gravito-turbulent model, however, we also require an esti- 
mate of the cooling time to estimate ( = f CO oifi-|Rafikov (2007} cal- 
culate the approximate cooling time for a convective and radiative 
disk (depending on whether or not it is convective and, if so, ac- 
counting for the rate-limiting of cooling by the disk photosphere). 
This gives 



Eel 



0"SB 



J^4 



f(r) 



(35) 



:2x 10 4 yr 



10 3 gcm~ 



T y 3 f(r) 
100K7 10 3 



where f(r) is a function (shown in Appendix [X} of the opacity 
which interpolates between the optically thin/thick, and convec- 
tive/radiative regimes. 

With these parameters calculated, for a given assumption 
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about what drives the turbulence - e.g. MRI, gravitoturbulence, 
convection, etc. - the compressive Mach number M c can be calcu- 
lated following §[6] We also technically need to assume the details 
of the turbulent spectral shape, for which we will assume a spectral 
index p = 2 and non-intermittent T — 0, as well as the gas equa- 
tion of state, for which we take 7 = 7/5, appropriate for molecular 
hydrogen. But these choices have small effects on our results, as 
shown in Fig.[T] 

7.2 The Characteristic Initial Fragment Mass 

In Fig. [3] we use this model to calculate the expected mass of a 
self-gravitating "fragment." Varying Eo, 1000, oe, and M*, we cal- 
culate the expected T and Q, assuming a constant accretion rate 
of M = 3 x 10 -7 MQyr~' at all radiiFlwhich dominates the disk 
temperature inside ~ 10 au. Given this, we calculate the mass spec- 
trum as Fig.[T| and define an average (M co iiapse} (the mass-weighted, 
spectrum integrated mass). Technically this depends on M c hence 
the turbulent driving mechanism, but the dependence is weak so we 
just assume M. c — 0.1 in all cases. 

In each case, (M co ii a p Se } ~ p~ Q 2 Mdi S k(< r*) as expected. Since 
p = (it E r\ ) /M» , this increases with disk surface density or mass, 
and also with increasing disk-to-stellar mass ratio. Recall r c ff iacc oc 
(MJ7 2 ) 1 / 4 oc r, 3/4 and Tiff,* oc r* l/2 . So modulo opacity correc- 
tions we expect (Af co iiapse) oc r^ 1 ' 5-1 *' at small radii < a few au, 
weakly increasing with r*; and (M co iiapse) oc ,.o 5 + 3 ( L5 - a ) at i ar g e 
r* , increasing more rapidly. 

This is only the initial self-gravitating, bound mass - it may 
easily evolve in time, as discussed in §[8] However, it is interesting 
that there is a broad range of masses possible, with Earth and super 
Earth-like masses more common at < 1 au and giant planet masses 
more common at > 10 au. 

7.3 Disk Temperatures at Which Direct Collapse Occurs 

Given a mass profile, then for a source of turbulence in §[6] we can 
translate the criteria for statistical stability - a probability /fj" <C 1 
of forming a fragment in a characteristic timescale ~Myr - into a 
range of midplane temperatures T m a. 

Fig.|4]shows this for a disk with Eo. 1000 = 1 and a = 1, around 
a solar-type star, as well as a disk with Eo. iooo = 10. Choosing 
ct = 1.5 gives a similar result but with the curve slopes system- 
atically shifted. Together this spans a range in disk mass at ~ 10 au 
of ~ 0.004- 0.07 M». We compare the expected Tmi(r», E), for 
accretion rates M = and M = 3x 10~ 7 Mq yr _1 . 

There is some T(r») below which Q < 1, so the disk will 
catastrophically fragment. But this is quite restrictive: Q < 1 only 
when E and r, are large (for Eo, iooo ~ 10 — 100, r, > 10 — 100 au 
for a ~ 1 — 1.5, respectively; roughly where Mdisk(< ''*) > 0.1M*). 

If the disk is convectively unstable, the resulting Mach num- 
bers lead to a large temperature range over which Q > 1, so the 
disk is classically stable, but P£L ~ 1; this is less likely to be rele- 
vant for a MMSN (r* > 50 au for Eo. iooo = 1) but can be sufficient 
at r» > 2 — 10 au for Eo, iooo > 10. As discussed in § |6.3| the con- 
vective Mach numbers are somewhat uncertain, so we show the 
calculation for the range therein. 

Strong-field MRI - if/where it is active - produces even larger 
Ml; this can greatly expand the range where _P/ r "' g ~ 1 even in a 
MMSN. Again there is a range of possible M c in the saturated 

7 This may not be self-consistent, since M could vary with disk parame- 
ters, but there is no straightforward a priori expectation for M, and we only 
intend this as a guide, in any case. 



state, shown here. In the weak-field regime, M c is smaller, giving 
results very similar to the convection prediction. 

Gravito-turbulence (again, if/where it is active), interestingly, 
has the opposite dependence on temperature. Because cooling rates 
grow rapidly at higher T^a, the expected Air is larger and P^' is 
larger at higher T (despite higher Q). It is also much less sensitive 
to the disk surface density. If the process operates, this is sufficient 
to produce ~ 1 at r* > 2 — 5 au, regardless of M for nearly all 
reasonable temperatures, and even at r* < 1 au if there is a modest 
M to raise the temperature (hence cooling rate) to > 300 — 1000 K. 
For the disk surface densities here, P^L <C 1 at r* < 1 au requires 
M < 3 x 10 _7 MQyr _1 . However, we caution that at sufficiently 
high T and Q, the mechanism may not operate at all. 

Beyond a certain radius (which depends on which combina- 
tion of these mechanisms are active), perhaps the most important 
result is that there may be no temperature where P"^ g <C 1, for a 
given surface density. 

7.4 Surface Densities at Which Direct Collapse Occurs 

In Fig. [5] we calculate the surface densities (as a function of radius 
r» around a solar-type star) where disks are statistically unstable 
(Pfng ~ 1) if/when different turbulent driving mechanisms are ac- 
tive, as in Fig. [4] Whereas in Fig.|4]we allowed the temperature to 
be free, here we assume it follows our best estimate T m a (r* , E) for 
either an accretion rateM = or M = 3 x 10 _7 Mq yr _1 , but freely 
vary E(r„). 

Roughly speaking, convection and MRI systematically lower 
the density at all radii where Pj£L ~ 1. Classical instability (Q < 1) 
requires E > 30Emmsn- If the disk is convective it can have Q > 1 
but be statistically unstable (vulnerable to direct fragmentation via 
turbulence density fluctuations) for E > IOEmmsn at r, > lOau; at 
smaller radii the threshold is sensitive to accretion-heating raising 
Q. If strong-field MRI is active, the threshold surface density for 
statistical instability is much lower: for small accretion rates and/or 
large radii > 2 — 5 au, even the MMSN can have P^L ~ 1 . Gravito- 
turbulence, if active, generates sufficient turbulence for statistical 
instability ~ 1) even at E ~ 0.1 Emmsn, at radii > 1 au re- 
gardless of accretion rate, and < 1 au for M>3x 10~ 7 Mq yr _ 1 . 

8 DISCUSSION & CONCLUSIONS 

Traditionally, a disk with Toomre Q > 1 is classically stable against 
gravitational collapse on all scales (modulo certain global gravita- 
tional instabilities). However, we show here that this is no longer 
true in a turbulent disk. Random turbulent density fluctuations can 
produce locally self-gravitating regions that will then collapse, even 
in Q > 1 disks. Formally, the probability of such an event is al- 
ways non-zero, so strictly speaking turbulent disks are never "com- 
pletely" stable, but can only be so statistically, if the probability of 
forming a self-gravitating fluctuation is small over the timescale of 
interest. Moreover, we can analytically predict the probability, as a 
function of total self-gravitating mass, of the formation of such a 
region per unit time in a disk (or disk annulus) with given proper- 
ties. 

We do this using the excursion-set formalism developed in Pa- 
per I-Paper II, which allows us to use the power spectra of turbu- 
lence to predict the statistical properties of turbulent density fluc- 
tuations. In previous papers, this has been applied to the structure 
of the ISM in galactic or molecular cloud disks; however, we show 
it is straightforward to extend to a Keplerian, proto-planetary disk. 
The most important difference between the case here and a galac- 
tic disk is that in the galactic case, turbulence is highly super-sonic 
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(A4 ~ 10 — 100), as opposed to sub/trans-sonic. And in galactic 
disks, cooling is rapid and the disk is globally self-gravitating (non- 
Keplerian), so systems almost always converge rapidly to Q w 1 
(see Hopkins et al. 2012); here, we expect a wider range of Q. And 
finally, proto-planetary disks are very long-lived relative to their 
local-dynamical times, so even quite rare events (with a probability 
of, say, 10~ 6 per dynamical time) may be expected over the disk 
lifetime. 

At Q rs 1, disks with A4 C > 0.1 are classically stable but sta- 
tistically unstable: they are likely, over the lifetime of the disk, 
to experience at least an order-unity number of "fragmentation" 
events (formation of self-gravitating, collapsing masses). As ex- 
pected, higher-<2 disks are "more stable" (although again we stress 
this is only a probabilistic statement); for Q ~ 3 — 5, values of 
M c > 0.3 — 0.5 are required for statistical instability. If the turbu- 
lence is transsonic (A4 C ~ 1 — 3), values as large as Q ~ 40 — 1000 
can be statistically unstable! Generally speaking, we show that sta- 
tistical stability (i.e. ensuring that the probability of a stochastic 
direct collapse event is <C 1) integrated over some timescale of in- 
terest to requires a Qmm 

2mm w 0.5 exp [ s/2 In (to JtrHl) In ( 1 + Ml)} 
~ 0.5 exp(6 v / MT+A*I)) 

(Fig. [2]& Eq. |16[>. At the Mach numbers of interest, this is expo- 
nentially increasing with M ! 

This is a radical revision to traditional stability criteria. How- 
ever, the traditional Toomre Q criterion is not irrelevant. It is a 
necessary, but not sufficient, criterion for statistical stability, which 
should not be surprising since its derivation assumes a homoge- 
nous, non-turbulent disk. If Q <C 1, then "catastrophic" fragmen- 
tation occurs even for M c — > 0; all mass in the disk is (classi- 
cally) unstable to self-gravity and collapse proceeds on a single 
free-fall time. If Q > 1, fragmentation transitions to the stochastic 
(and slower) statistical mode calculated here, dependent on ran- 
dom turbulent density fluctuations forming locally self-gravitating 
regions. 

Likewise, the criterion injGammie (2 001} , that the cooling 
time be longer than a couple times the dynamical time, is not suf- 
ficient for statistical stability. For a given turbulent Mach number 
and g, we show that assuming a stiffer equation of state has quite 
weak effects on the "stochastic" mode of fragmentation; in fact, 
even pure adiabatic gas (no cooling) produces very similar statis- 
tics if it can sustain a similar M c - Again, the key is that the jGam-] 
|mie| ( |200T| > criterion is really about the prevention of catastrophic 
fragmentation; as noted therein, a sufficiently slow cooling time al- 
lows turbulence to maintain a steady state Q ~ 1, and dissipation of 
that turbulence (driven by gravitational density waves and inflow) 
can maintain the gas thermal energy (c,). Thus faster cooling leads 
to catastrophic fragmentation of most of the mass on a single dy- 
namical time (Q < 1 and M c 3> 1). And indeed this is the case 
from Paper I in a galactic disk, where f coo i <C /d yn and the mass is 
only "recycled" back into the diffuse medium by additional energy 
input (from stellar feedback). 

However, provided the Gammie (2001) criterion is met and 
Q > 1 everywhere, we still predict fragmentation in the "stochas- 
tic" mode if gravito-turbulence operates. The cooling time is then 
important insofar as it changes the equilibrium balance between 
turbulent and thermal energy, i.e. appears in Q and governs M c oc 
('cool f2)~ 1//2 . This, indeed, has now been seen in a growing number 
of simulations (see references in § [TJ, involving either larger vol- 
umes and/or longer runtimes. We consider the application of our 



theory to these specific models in § |6. 1| this allows us to predict a 
revised cooling-time criterion in this "mode," required for statisti- 
cal stability over any timescale of interest: 

c = jcooi^ 4 ln[f ^- 2 O)] 

> 3 7 (7- 1) [ln(2Q)P ' 

This provides an excellent explanation for the results in these sim- 
ulations, and resolves the apparent discrepancies between them 
noted in §[T] 

If another process is able to drive turbulence, then stochas- 
tic direct collapse might occur with even longer cooling times. We 
show that if the MRI is active and saturates at strong-field /3 sat ~ 1, 
the required Q for complete suppression of fragmentation can be 
very large (> 10— 15; scaling with /3 silt as Eq.|28|), independent of 
the cooling time. If the MRI is not active (in the dead zone, for ex- 
ample), or if it saturates at weak-field values /3 sat > 10, and cooling 
is slower than the limit above, then convection may be the dom- 
inant driver of turbulence. This is sufficient to produce stochastic 
fragmentation events in the range Q ~ 1 — 3, though probably not 
much larger. 

We apply these calculations to specific models of proto- 
planetary disks that attempt to self-consistently calculate their tem- 
peratures and cooling rates. Doing so, we show that the parame- 
ter space where statistical instability and stochastic fragmentation 
may occur is far larger than that of classical instability (where 
Q < 1), and can include most of the disk even in a MMSN. 
Gravito-turbulence appears to be the most important channel driv- 
ing stochastic fragmentation when it is active, and is sufficient 
to produce an order-unity number of events in disks with E > 
0. 1 Emmsn at distances > 1 au (independent of accretion rate) and 
even at < 1 au (if the disk is heated by accretion rates > 3 x 
10 -7 MQyr~'). At low accretion rates and/or large radii, strong- 
field MRI (if active) is also sufficient to drive statistical instability 
if E > Emmsn- And we show that beyond a few au, the combination 
of gravito-turbulence, convection, and Toomre instability at low-T 
means there may be no disk temperatures at which f& < 1 in a 
modest-density disk. 

Ultimately, regardless of our (admittedly speculative) discus- 
sion of theoretical models for the sources of turbulence in proto- 
planetary disks, the key question is empirical. Are the actual Mach 
numbers in such disks sufficient, for their Q values, to be "inter- 
esting" here? This remains an open question. However, Hughes 
|et al.| ( [201 1\ present some early indications of detection of tur- 
bulent linewidths in two protoplanetary disks, with inferred Mach 
numbers of ~ 0.1 and ~ 0.4. Although uncertain, these essen- 
tially bracket the most interesting regime of our calculations here! 
Because of the exponential dependence of stochastic collapse on 
Mach numbers, future observations which include larger statisti- 
cal samples and more/less massive disks, as well as constraints on 
whether the turbulence appears throughout the disk (since it is the 
midplane Mach numbers that matter most for the models here), will 
be critical to assess whether the processes described in this paper 
are expected to be relatively commonplace or extremely rare. 

We also predict the characteristic mass spectrum of fragmen- 
tation events. Sub-sonic turbulence produces a narrow mass spec- 
trum concentrated around the Toomre mass ~ /1 2 Q 2 Ma^; angular 
momentum and shear suppress larger-scale collapse while thermal 
(and magnetic) pressure suppress the formation of smaller-scale 
density fluctuations. For typical mass ratios fi ~ 0.1 in the early 
stages of disk evolution, this corresponds reasonably well to the 
masses of giant planets. For smaller mass ratios fi ~ 0.01, which 
should occur at somewhat later stages of evolution, this implies 
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direct collapse to Earth-like planet masses may be possible. Of 
course, such collapse will carry whatever material is mixed in the 
disk (i.e. light elements), so such a planet would presumably lose its 
hydrogen/helium atmosphere as it subsequently evolved (see refer- 
ences in §[TJ. As the turbulence approaches transsonic, the mass 
spectrum becomes much more broad. As shown in Paper II, this 
owes to the greater dynamic range in which turbulence is impor- 
tant; in the limit A4 — > oo, the mass spectrum approaches a power 
law with equal mass at all logarithmic intervals in mass (up to the 
maximum disk Jeans mass). Thus within a more turbulent disk, di- 
rect collapse to a wide range of masses - even at identical disk 
conditions and radii - is expected. This may explain observed sys- 
tems with a range of planet masses within narrow radii (e.g. [Carter] 
|et al.|2012 ), and it also predicts a general trend of increasing av- 
erage collapse mass with distance from the central star. However, 
we are not accounting for subsequent orbital evolution here, and 
subsequent accretion (after collapse) will modify the mass spectra. 

We show how these predictions can be modified to allow for 
different turbulent velocity spectra, equations-of-state of the gas, 
the presence/absence of magnetic fields, and different disk mass 
profiles; however these are often small corrections or amount to 
order-unity re-normalizations of the predicted mass spectra. They 
do not change any of our qualitative conclusions. 

Throughout, we restrict our focus to the formation of self- 
gravitating regions. Following the subsequent evolution of those 
regions (collapse, fragmentation, migration, accretion, and possi- 
ble formation of planets) requires numerical simulations to treat 
the full non-linear evolution (see e.g. Kratter & Murray-Clay 201 1 



|Rice et al.|[20TT| |Zhu et alpOH) |Galvagni et al.||2012| and ref- 
erences therein). Ideally, this would be within global models that 
can self-consistently follow the formation of these regions. How- 
ever, this is computationally extremely demanding. Even ignoring 
the detailed physics involved, the turbulent cascade must be prop- 
erly followed (always challenging), and most important, since the 
fluctuations of interest can be extremely rare, very large (but still 
high-resolution) boxes must be simulated for many, many dynami- 
cal times. As discussed in §[T] this has led to debate about whether 
or not different simulations have (or even can be) converged. Most 
of the longest-duration simulations to date have been run for times 
~ lOOOfi - , which for plausible M c and Q may be shorter by a 
factor of ~ 10 3 — 10 5 than the timescale on which of order a single 
event is expected to occur in the entire disk. But certainly in the 
example of planet formation, a couple of rare events are "all that is 
needed," so this is an extremely interesting case for future study. 
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Figure Al. Characteristic mass at collapse, as Fig. [5] but with a more de- 
tailed set of opacity tables and temperature calculation as described in § [A| 



APPENDIX A: DETAILS OF TEMPERATURE AND 
COOLING RATE CALCULATION 

Recall, from the text, in the case of disk irradiation by a central 
solar-type star, the effective temperature is 

r eff ,«^V /4 r* 



1997 1 



I Chiang & Goldreich 
0.05 (n/smY'^R, 



this gives an approximate r m id. 



(Al) 

0.005 (r„/au)-' + 



with a 

Rq and T» = 6000 A". In this regime the ex- 
ternal radiation produces a hot surface dust layer which re-radiates 
~ 1/2 the absorbed light back into the disk, maintaining 7^; if 
the disk is optically thick to the incident and re-radiated emission, 

C/2' /4 



Meanwhile, accre- 



tion produces an effective temperature 

3 Mf2 z 



Toft . a 



r 3 MWi 

L87T (JSB J 



1/4 



(A2) 



In the optically thick limit, this is just related to T m \& by T* M acc ~ 
(3/4) {tr + 2/3) r c f f _ acc , where tr is the Rosseland-mean optical 
depth, t r — KR(T m i<i) E (so r m id. aC c is determined implicitly). 

In the text, Figs. |4|5| simply use the optically thick relations 
for Zlnid and, when accretion is included, interpolate with T^ ld ~ 
?mid * + Tmid. acc • We approximate the opacities with the simple k r ~ 
5cm 2 g"' at T > 160 K, and k r ~ 2.4 x 10~ 4 T 2 cnrg~' K~ 2 at 
lower temperatures (approximately what is obtained with a galactic 
gas-to-dust ratio; se e|Adams & Shu|1988||Bell & Lin|1994) . 

We then follow Rafikov j2007) andestimate the cooling time 

as 



^cool 



f(r) 



(A3) 



where c, = yj kg T m \i / fi with fi ~ 2.3 appropriate for molecular, 
dusty gas. The function fir) is given by the interpolation between 
the convective and radiative terms at the photosphere 

-i 



f{r) 



ry- 



■ XT v + ( 



4 7 - 1 ( 7 -l) 



(A4) 
(A5) 



l+a + /S7- 1 (7-l) 

here \ and <f> are constants; from the detailed estimates therein 
(f> ~ 1 and x ~ 0.19 — 0.31, depending on the disk parameters, 
so we assume \ = 0.31 to be conservative (since this gives larger 
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Figure A2. Shaded regions show the temperature range in which there is 
statistical instability (order-unity probability of a collapse event), as Fig. [4] 
but with a more detailed set of opacity tables and temperature calculation 
as described in §[A] 



cooling times). This interpolates between the optically thick lim- 
its (dominated by x tV ) an d optically thin cases (4>t~\ where the 
cooling flux becomes oc tosb T^ id ). 

The scaling index r\ is determined from the temperature gradi- 
ent to the photosphere for a disk in vertical hydrostatic equilibrium, 
under the assumption that over some (limited) local temperature 
range the opacity k can be approximated by k ~ reo P a T 13 , which is 
valid for our assumptions. We assume 7 = 7/5 in both this and the 
turbulent density fluctuation calculation in Figs. |3|5| Based on the 
scaling of opacities in Sem enov et al.|p003) , at T < 160 K, a = 0, 
(J«2, so 77 » 7/11; at 160 < T < 1500K a = 0, /3 » 1/2- 1 and 
77 ~ 7/8 — 7/9 (we adopt 7/9, but this makes no significant differ- 
ence), and at T > 1500 K (when all grains sublimate and molecular 
opacity dominates) a ~ 2/3 and /3 « 7/3, so r] ~ 3/7. 

In the text, we restrict to these simple estimates because the 
quantities of interest are fairly uncertain. We can, however, examine 
a more detailed approximation here. First, we take the opacities 
/t«(?mid) from the full tabulated values in Semenov et al. 1 2003 1. 



Second, a more accurate estimate of the midplane temperature is 
given by solving the implicit equation 



T 4 - 



2 

3tv J 



~» + - + - — ] Tau 



cc+ 1 + 



(A6) 



where tv = tv (T m n) is the vertical optical depth from the midplane, 
tv = K«(7mid) S/2. This allows for an appropriate interpolation be- 
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Figure A3. Shaded regions show the surface density range in which there is 
statistical instability (order-unity probability of a collapse event), as Fig. [5] 
but with a more detailed set of opacity tables and temperature calculation 
as described in §|Al 



tween the optically thick case and the case where disk is optically 
thin to its own re-radiation. Third, we can switch between the /(r) 
above, appropriate for a convective disk, and /(r) w r + t~' ap- 
propriate for a purely radiative, convectively stable disk, when the 
disk falls below the (temperature-dependent) criteria for convec- 
tive instability V > V a d with V = (1 + q)/(4 - /3 ) and V ad = 
(7 — 1 )/7 (see Lin & Papaloizou 1980 Rafikov 2007 1, where a and 



/3 depend on T (using the full explicit derivatives from the opacity 
tables). 

Figures fA 1|A3 | repeat our calculations from the main text with 
this more detailed temperature calculation. We find that the quanti- 
tative results in Figs. |3|5| are all changed at the factor < 2 level and 
the qualitative conclusions are completely unchanged. The sense of 
the quantitative change tends to slightly expand the regions of pa- 
rameter space where statistical instability and turbulence-promoted 
fragmentation can occur. The more detailed opacity calculation im- 
prints some small features on the parameter space, the most signif- 
icant of which is the elimination of predicted temperatures much 
larger than ~ 1500 K at small radii (because grains sublimate and 
cooling becomes optically thin), but the "gravito-turbulent" thresh- 
olds shift with the predicted temperatures so the statistical stability 
is essentially identical. 
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